library(tidyverse)
library(here)
library(haven)
library(fastDummies)
library(glmnet)
library(modelsummary)
library(ggpubr)
library(estimatr)
library(wesanderson)
library(Matrix)

rm(list=ls())

## program to run lasso with rhs and dv provided. 
get_nonzero_coefficients <- function(rhs, dv) {
  lasso <- glmnet(rhs, dv)
  lambda_1se <- lasso$lambda[which(lasso$lambda >= mean(lasso$lambda) - sd(lasso$lambda))]
  lambda_1se <- lambda_1se[which.min(abs(lambda_1se - mean(lasso$lambda)))]
  coefmat <- as.matrix(coef(lasso, s = lambda_1se))
  coefmat <- as.data.frame(coefmat)
  nonzero <- coefmat %>% filter(s1 > 0)
  row_names_vector <- as.vector(row.names(nonzero))
  row_names_vector <- row_names_vector[-1]
  return(row_names_vector)
}

## Load cleaned and anonymous data:
d<-readRDS(here("Data","cleandata.rds"))
  
  
## Tidy up the data a bit more: 
d <- d %>%  mutate(petition_ordinal = recode(petition, "Yes" = 2, "Maybe" = 1, "No" = 0))


### create matrix of control variables for the lasso selection: 
vars<- d %>% select(pool_concede, pool_concede_notstd, treat, info, certain, uncertain, control,
                    debtpreA, debtpreB, debtpreC,
                    partyid, age, Q9, inc, education,
                    approval_biden, RWnews, urbrur, passcheck,
                    wait_pool, blame_mult, signpet,ceiling_import)
vars <- vars %>% drop_na(pool_concede)

vars <- vars %>%
  mutate(education = recode(education,
                            "Associates or technical degree" = "some college",
                            "Bachelor’s degree" = "college graduate",
                            "Graduate or professional degree (MA, MS, MBA, PhD, JD, MD, DDS etc.)" = "college graduate",
                            "High school diploma or GED" = "high school or less",
                            "Prefer not to say" = "high school or less",
                            "Some college, but no degree" = "some college",
                            "Some high school or less" = "high school or less"))

vars <- vars %>%
  mutate(education = recode(education,
                            "high school or less" = 1,
                            "some college" = 2,
                            "college graduate" = 3))

## replace missing covariates with mean or median 
vars <- vars %>%
  mutate(inc = ifelse(is.na(inc), median(inc, na.rm = TRUE), inc))

vars <- vars %>%
  mutate(urbrur = ifelse(is.na(urbrur), median(urbrur, na.rm = TRUE), urbrur))

vars <- vars %>%
  mutate(approval_biden = ifelse(is.na(approval_biden), median(approval_biden, na.rm = TRUE), approval_biden))

vars <- vars %>%
  mutate(female_dummy = ifelse(Q9 == "Female", 1, 0))

rhs <- vars %>% select(debtpreA,debtpreB,debtpreC,partyid,age,female_dummy,inc,education,approval_biden,RWnews,urbrur,passcheck)

## Create dummies to include in the LASSO model
rhs <- rhs %>% fastDummies::dummy_cols("partyid")
rhs <- rhs %>% fastDummies::dummy_cols("inc")
rhs <- rhs %>% fastDummies::dummy_cols("education")
rhs <- rhs %>% fastDummies::dummy_cols("RWnews")
rhs <- rhs %>% fastDummies::dummy_cols("approval_biden")
rhs <- rhs %>% fastDummies::dummy_cols("urbrur")

sumvars <- as.data.frame(rhs)
dvsum <- vars %>% select(pool_concede_notstd)
treatsum <- vars %>% select(certain, uncertain, info)
sumvars <- cbind(sumvars,dvsum,treatsum)
rm(dvsum, treatsum)

rhs <- rhs %>%
  select(-partyid, -inc, -education, -RWnews, -approval_biden, -urbrur )

rhs<-as.matrix(rhs)
dv <- as.matrix(vars$pool_concede)

row_names_vector <- get_nonzero_coefficients(rhs, dv)

vars <- cbind(vars,rhs)


#### Main Models: 


### Comparison of Info vs. No Info 
formula <- as.formula(paste("pool_concede ~ info + ", paste(row_names_vector, collapse = "+")))
lm_info <- lm_robust(formula, d=vars)
modelsummary(lm_info, stars=TRUE)

# Comparison of Certain and Uncertain vs. No Info 
formula <- as.formula(paste("pool_concede ~ certain + uncertain +", paste(row_names_vector, collapse = "+")))
lm_certain <- lm_robust(formula, d=vars)
modelsummary(lm_certain, stars=TRUE)

# Comparison of Certain vs. Uncertain 
formula <- as.formula(paste("pool_concede ~ certain + control +", paste(row_names_vector, collapse = "+")))
lm_certain2 <- lm_robust(formula, d=vars)
modelsummary(lm_certain2, stars=TRUE)


# Table for Supplementary Appendix 
models<-list(
  "Information vs. Control" = lm_info, 
  "Certain & Uncertain vs. Control" = lm_certain, 
  "Certain vs. Uncertain" = lm_certain2
)

cmap <- c("info" = "Information", "certain"="Certain", "uncertain"="Uncertain", "control"="Control",
          "debtpreA" = "Debt General",
          "female_dummy" = "Woman",
          "partyid_2" = "Party ID (2)",
          "inc_4" = "Income (4)",
          "education_3" = "Education (3)",
          "approval_biden_2" = "Biden Approval (2)")

notelist <- paste(row_names_vector, collapse = ", ")
note <- paste("Non-zero covariates include: ", notelist, collapse="")
modelsummary(models, stars=T, notes=note, coef_map=cmap)
#modelsummary(models, stars=T, notes=note , output="latex", coef_map=cmap)



### Balance Table: 
sumvars <- sumvars %>% select(debtpreA, debtpreB, debtpreC ,
                              partyid,age,female_dummy,inc,
                              education,approval_biden,RWnews,
                              urbrur,info,certain,uncertain)
sumvars <- as.data.frame(sumvars)

sumvars <- sumvars %>% rename(`Debt General` = debtpreA,
                              `Debt Spending` = debtpreB,
                               `Debt Taxes` = debtpreC,
                              PartyID = partyid,
                              Age = age,
                              Woman = female_dummy,
                              Income = inc,
                              Education = education,
                              `Biden Approval` = approval_biden,
                              `Right Wing News` = RWnews,
                              `Urban Rural` = urbrur
                              )

datasummary_skim(sumvars, type="numeric")
datasummary_balance(~info, data = sumvars)
#datasummary_balance(~info, data = sumvars, output = 'latex')





#### simple histogram from ANOVA. 

dvtreat <- vars %>% select(pool_concede_notstd,treat)
dvtreat$treat2 <- as.factor(dvtreat$treat)
#### "simple" histogram from ANOVA.

means<- dvtreat %>% 
  group_by(treat2) %>% 
  summarise(
    mean_sleep = mean(pool_concede_notstd),
    sd_sleep   = sd(pool_concede_notstd)  ) 


means <- means %>%
  mutate(treat2 = recode(treat2, "0" = "Control", "1" = "Certain", "2" = "Uncertain"))

p_value_one <- tibble(
  x = c("Control", "Control", "Certain", "Certain"),
  y = c(3.5, 3.6, 3.6, 3.5)
)

p_value_two <- tibble(
  x = c("Control", "Control", "Uncertain", "Uncertain"),
  y = c(3.0, 3.1, 3.1, 3.0)
)

p_value_three <- tibble(
  x = c("Certain", "Certain", "Uncertain", "Uncertain"),
  y = c(2.5, 2.6, 2.6, 2.5)
)



#c01 <- lm_robust(pool_concede_notstd~treat2, data=dvtreat)
modelsummary(lm_certain, stars=T)
p1<- round(lm_certain$p.value["certain"], 7)
p1 <- ifelse(p1<.0001,0.000,p1)
p1 <-pvalue_string <- paste("p=", p1)
p1 <- ifelse(p1=="p= 0","p= 0.000",p1)
p1

#fx size
lm_certain$coefficients["certain"]/sd(dvtreat$pool_concede_notstd)
lm_certain$coefficients["uncertain"]/sd(dvtreat$pool_concede_notstd)

p2 <- lm_certain$p.value["uncertain"]
p2<- round(p2, 3)
p2 <- ifelse(p2<.0005,0.000,p2)
p2 <-pvalue_string <- paste("p=", p2)
p2 <- ifelse(p1=="p= 0","p= 0.000",p2)
p2



p3<- round(lm_certain2$p.value["certain"], 3)
p3 <- ifelse(p3<.0005,0.000,p3)
p3 <-pvalue_string <- paste("p=", p3)
p3 <- ifelse(p1=="p= 0","p= 0.000",p3)
p3

barplot <- means %>% 
  ggplot(aes(treat2, mean_sleep)) +
  geom_col(aes(fill = treat2), color = "black", width = 0.85) +
  geom_line(data = p_value_one, 
            aes(x = x, y = y, group = 1)) +
  geom_line(data = p_value_two, 
            aes(x = x, y = y, group = 1)) +
  geom_line(data = p_value_three, 
            aes(x = x, y = y, group = 1)) +
  annotate("text", x = 1.5, y = 3.71, 
           label = "p=0.000",
           size = 3, color = "#22292F") +
  annotate("text", x = 2, y = 3.21, 
           label = p2,
           size = 3, color = "#22292F") +
  annotate("text", x = 2.5, y = 2.71, 
           label = p3,
           size = 3, color = "#22292F") +
  scale_fill_manual(values = wes_palette("Royal1")) +
scale_y_continuous(limits = c(0, 4), expand = c(0, 0)) +
  guides(fill = FALSE) +
  ggtitle("(b) ") + theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Treatment Condition", y = "") +
  theme(panel.background = element_rect(fill = "grey94"))


barplot 

#### "simple" histogram from ANOVA. binary comparison 
vars$info2 <- as.factor(vars$info)
dvtreat <- vars %>% select(pool_concede_notstd,info2)

means<- dvtreat %>% 
  group_by(info2) %>% 
  summarise(
    mean_sleep = mean(pool_concede_notstd),
    sd_sleep   = sd(pool_concede_notstd)  ) 

means <- means %>%
  mutate(info2 = recode(info2, "0" = "Control", "1" = "Information"))


p_value_one <- tibble(
  x = c("Control", "Control", "Information", "Information"),
  y = c(2.5, 2.6, 2.6, 2.5)
)

modelsummary(lm_info, stars=T)
p1<- round(lm_info$p.value["info"], 3)

p1 <- ifelse(p1<.0001,0.000,p1)
p1 <-pvalue_string <- paste("p=", p1)
p1 <- ifelse(p1=="p= 0","p= 0.000",p1)


#fx size: 
lm_info$coefficients["info"]/sd(dvtreat$pool_concede_notstd)

barplot_2 <- means %>% 
  ggplot(aes(info2, mean_sleep)) +
  geom_col(aes(fill = info2), color = "black", width = 0.85) +
  geom_line(data = p_value_one, 
            aes(x = x, y = y, group = 1)) +
  annotate("text", x = 1.5, y = 2.71, 
           label = "p=0.000",
           size = 3, color = "#22292F") +

  scale_fill_manual(values = wes_palette("Royal1")) +
  scale_y_continuous(limits = c(0, 4), expand = c(0, 0)) +
  guides(fill = FALSE) +
  ggtitle("(a)") + theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Treatment Condition", y = "") + 
  theme(panel.background = element_rect(fill = "grey94"))
barplot_2



bplots <- ggarrange(barplot_2,barplot,label.x = "Treatment Conditions")
bplots
#ggsave("pool_barplot.jpg",bplots)


### Comparison of Info vs. No Info 
formula <- as.formula(paste("pool_concede ~ info + ", paste(row_names_vector, collapse = "+")))
lm_info <- lm_robust(formula, d=vars)
modelsummary(lm_info, stars=TRUE)

# Comparison of Certain and Uncertain vs. No Info 
formula <- as.formula(paste("pool_concede ~ certain + uncertain +", paste(row_names_vector, collapse = "+")))
lm_certain <- lm_robust(formula, d=vars)
modelsummary(lm_certain, stars=TRUE)

vars$certain2 <- vars$certain

# Comparison of Certain vs. Uncertain 
formula <- as.formula(paste("pool_concede ~ certain2 + control +", paste(row_names_vector, collapse = "+")))
lm_certain2 <- lm_robust(formula, d=vars)
modelsummary(lm_certain2, stars=TRUE)

models<-list(
  "Information vs. Control" = lm_info, 
  "Certain & Uncertain vs. Control" = lm_certain, 
  "Certain vs. Uncertain" = lm_certain2
)

cmap <- c('certain2' = 'Certain (base Uncertain)', 'certain' = 'Certain (base Control)', 'uncertain' = 'Uncertain (base Control)', 'info' = 'Information (base Control)')

coefplot <- modelplot(models, coef_map = cmap) + labs(x = " ")

coefplot2 <- coefplot + 
  scale_color_manual(values = c("#C93312", "#C93312", "#C93312", "#C93312")) + 
  theme(panel.background = element_rect(fill = "grey95"),
        panel.grid.major = element_line(color = "white", size = 0.5),
        panel.grid.minor = element_line(color = "white", size = 0.25)) + 
  ggtitle("(c)") + 
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "none",  # Remove the legend
        legend.background = element_rect(fill = "white", colour = "black"), 
        legend.key = element_rect(fill = "white", colour = "white")) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black")

coefplot2
plots3 <- ggarrange(bplots,coefplot2, ncol = 1)
plots3
#ggsave("pool_barplot_coefplot.jpg",plots3)




#### #### #### #### #### #### #### #### #### 
#### Dems and Republicans: 
#### #### #### #### #### #### #### #### #### 

vars<- d %>% select(pool_concede, pool_concede_notstd, treat, info, certain, uncertain, control,
                    debtpreA, debtpreB, debtpreC,
                    partyid, age, Q9, inc, education,
                    approval_biden, RWnews, urbrur, passcheck,
                    wait_pool, blame_mult, signpet)

vars <- vars %>% drop_na(pool_concede)
summary(vars)

dvtreat <- vars %>% select(pool_concede_notstd,treat,partyid,certain,uncertain,control)
dvtreat$treat2 <- as.factor(dvtreat$treat)
#### "simple" histogram from ANOVA.

dvtreat <- dvtreat %>% 
  mutate(dem = ifelse(partyid %in% c(0, 1), 1, 0))

dvtreat <- dvtreat %>% 
  mutate(rep = ifelse(partyid %in% c(5, 6), 1, 0))

dvtreat <- dvtreat %>% 
  mutate(ind = ifelse(partyid %in% c(2, 3, 4), 1, 0))


means<- dvtreat %>% 
  group_by(treat2,dem,rep,ind) %>% 
  summarise(
    means = mean(pool_concede_notstd),
    sd_dem   = sd(pool_concede_notstd),
    count = sum(!is.na(pool_concede_notstd)))

demrepcount <- dvtreat %>% 
  group_by(dem,rep,ind) %>% 
  summarise(
    count = sum(!is.na(pool_concede_notstd)))
    
means <- means %>%
  mutate(treat2 = recode(treat2, "0" = "Control", "1" = "Certain", "2" = "Uncertain"))

means <- means %>%
  mutate(repdem = case_when(
    rep == 1 ~ "Republican",
    dem == 1 ~ "Democrat",
    ind == 1 ~ "Independent",
    TRUE ~ "Other"
  ))

p_value_one <- tibble(
  x = c("Control", "Control", "Certain", "Certain"),
  y = c(3.5, 3.6, 3.6, 3.5))

p_value_two <- tibble(
  x = c("Control", "Control", "Uncertain", "Uncertain"),
  y = c(3.1, 3.2, 3.2, 3.1))

p_value_three <- tibble(
  x = c("Certain", "Certain", "Uncertain", "Uncertain"),
  y = c(2.7, 2.8, 2.8, 2.7))




#dem_means<- means %>% filter(repdem=="Democrat")

temp<-cbind(rhs,dv)
temp<- as.data.frame(temp)

temp_dem <- temp %>% filter(partyid_0==1 | partyid_1==1)
temp_dem <- temp_dem %>% rename(pool_concede = V31)

dv_dem <- as.matrix(temp_dem$pool_concede)
rhs_dem <- temp_dem %>% select(-pool_concede)
rhs_dem <- as.matrix(rhs_dem)


row_names_vector <- get_nonzero_coefficients(rhs_dem, dv_dem)


vars_dem <- vars %>% filter(partyid<2)
vars_dem <- cbind(vars_dem,rhs_dem)
### Comparison of Info vs. No Info 
formula <- as.formula(paste("pool_concede_notstd ~  info +", paste(row_names_vector, collapse = "+")))
dem0 <- lm_robust(formula, d=vars_dem)
modelsummary(dem0, stars=TRUE)

formula <- as.formula(paste("pool_concede_notstd ~  certain + uncertain +", paste(row_names_vector, collapse = "+")))
dem1 <- lm_robust(formula, d=vars_dem)
modelsummary(dem1, stars=TRUE)


formula <- as.formula(paste("pool_concede_notstd ~  certain + control +", paste(row_names_vector, collapse = "+")))
dem2 <- lm_robust(formula, d=vars_dem)
modelsummary(dem2, stars=TRUE)

demnotes <- paste(row_names_vector, collapse = ", ")


#c01 <- lm_robust(pool_concede_notstd~treat2,data=dvtreat_dem )
#modelsummary(c01, stars=T)
p1<- round(dem1$p.value["certain"], 7)
p1 <- ifelse(p1<.0001,0.000,p1)
p1 <-pvalue_string <- paste("p=", p1)
p1 <- ifelse(p1=="p= 0","p= 0.000",p1)

#fxsize 
dem1$coefficients["certain"]/sd(vars_dem$pool_concede)
dem1$coefficients["uncertain"]/sd(vars_dem$pool_concede)

p2<- round(dem1$p.value["uncertain"], 4)
p2 <- ifelse(p2<.0001,0.000,p2)
p2 <-pvalue_string <- paste("p=", p2)
p2 <- ifelse(p2=="p= 0","p= 0.000",p2)



p3<- round(dem2$p.value["certain"], 4)
p3 <- ifelse(p3<.0001,0.000,p3)
p3 <-pvalue_string <- paste("p=", p3)
p3 <- ifelse(p3=="p= 0","p= 0.000",p3)


dem_barplot <- means %>% filter(repdem=="Democrat") %>%
  ggplot(aes(treat2, means)) +
  geom_col(aes(fill = treat2), color = "black", width = 0.85) +
  #  geom_errorbar(aes(ymin = mean_sleep - sd_sleep,
  #                    ymax = mean_sleep + sd_sleep),
  #                color = "#22292F",
  #                width = .1) +
  geom_line(data = p_value_one, 
            aes(x = x, y = y, group = 1)) +
  geom_line(data = p_value_two, 
            aes(x = x, y = y, group = 1)) +
  geom_line(data = p_value_three, 
            aes(x = x, y = y, group = 1)) +
  annotate("text", x = 1.5, y = 3.81, 
           label = "p=0.000",
           size = 3, color = "#22292F") +
  annotate("text", x = 2, y = 3.31, 
           label = p2,
           size = 3, color = "#22292F") +
  annotate("text", x = 2.5, y = 2.91, 
           label = p3,
           size = 3, color = "#22292F") +
  scale_fill_manual(values = wes_palette("Royal1")) +
  scale_y_continuous(limits = c(0, 4), expand = c(0, 0)) +
  guides(fill = FALSE) +
  ggtitle("Democrats") + theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Treatment Condition", y = "")  # Add this line to set the x-axis label

dem_barplot_notitle <- means %>% filter(repdem=="Democrat") %>%
  ggplot(aes(treat2, means)) +
  geom_col(aes(fill = treat2), color = "black", width = 0.85) +
  geom_line(data = p_value_one, 
            aes(x = x, y = y, group = 1)) +
  geom_line(data = p_value_two, 
            aes(x = x, y = y, group = 1)) +
  geom_line(data = p_value_three, 
            aes(x = x, y = y, group = 1)) +
  annotate("text", x = 1.5, y = 3.81, 
           label = "p=0.000",
           size = 3, color = "#22292F") +
  annotate("text", x = 2, y = 3.31, 
           label = p2,
           size = 3, color = "#22292F") +
  annotate("text", x = 2.5, y = 2.91, 
           label = p3,
           size = 3, color = "#22292F") +
  scale_fill_manual(values = wes_palette("Royal1")) +
  scale_y_continuous(limits = c(0, 4), expand = c(0, 0)) +
  guides(fill = FALSE) +
  ggtitle("(a)") + theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Treatment Condition", y = "")  # Add this line to set the x-axis label


dem_barplot 
dem_barplot_notitle


rep_means<- means %>% filter(repdem=="Democrat")

temp_rep <- temp %>% filter(partyid_5==1 | partyid_6==1)
temp_rep <- temp_rep %>% rename(pool_concede = V31)

dv_rep <- as.matrix(temp_rep$pool_concede)
rhs_rep <- temp_rep %>% select(-pool_concede)
rhs_rep <- as.matrix(rhs_rep)


row_names_vector <- get_nonzero_coefficients(rhs_rep, dv_rep)


vars_rep <- vars %>% filter(partyid>4)
vars_rep <- cbind(vars_rep,rhs_rep)
### Comparison of Info vs. No Info 
formula <- as.formula(paste("pool_concede_notstd ~  info +", paste(row_names_vector, collapse = "+")))
rep0 <- lm_robust(formula, d=vars_rep)
modelsummary(rep0, stars=TRUE)


formula <- as.formula(paste("pool_concede_notstd ~  certain + uncertain +", paste(row_names_vector, collapse = "+")))
rep1 <- lm_robust(formula, d=vars_rep)
modelsummary(rep1, stars=TRUE)

formula <- as.formula(paste("pool_concede_notstd ~  certain + control +", paste(row_names_vector, collapse = "+")))
rep2 <- lm_robust(formula, d=vars_rep)
modelsummary(rep2, stars=TRUE)

repnotes <- paste(row_names_vector, collapse = ", ")


p1<- round(rep1$p.value["certain"], 3)
p1 <- ifelse(p1<.0001,0.000,p1)
p1 <-pvalue_string <- paste("p=", p1)
p1 <- ifelse(p1=="p= 0","p= 0.000",p1)

p2<- round(rep1$p.value["uncertain"], 3)
p2 <- ifelse(p2<.0001,0.000,p2)
p2 <-pvalue_string <- paste("p=", p2)
p2 <- ifelse(p2=="p= 0","p= 0.000",p2)


#fxsize 
rep1$coefficients["certain"]/sd(vars_rep$pool_concede)
rep1$coefficients["uncertain"]/sd(vars_rep$pool_concede)



p3<- round(rep2$p.value["certain"], 3)
p3 <- ifelse(p3<.0001,0.000,p3)
p3 <-pvalue_string <- paste("p=", p3)
p3 <- ifelse(p3=="p= 0","p= 0.000",p3)


rep_barplot <- means %>% filter(repdem=="Republican") %>%
  ggplot(aes(treat2, means)) +
  geom_col(aes(fill = treat2), color = "black", width = 0.85) +
  #  geom_errorbar(aes(ymin = mean_sleep - sd_sleep,
  #                    ymax = mean_sleep + sd_sleep),
  #                color = "#22292F",
  #                width = .1) +
  geom_line(data = p_value_one, 
            aes(x = x, y = y, group = 1)) +
  geom_line(data = p_value_two, 
            aes(x = x, y = y, group = 1)) +
  geom_line(data = p_value_three, 
            aes(x = x, y = y, group = 1)) +
  annotate("text", x = 1.5, y = 3.81, 
           label = p1,
           size = 3, color = "#22292F") +
  annotate("text", x = 2, y = 3.31, 
           label = p2,
           size = 3, color = "#22292F") +
  annotate("text", x = 2.5, y = 2.91, 
           label = p3,
           size = 3, color = "#22292F") +
  scale_fill_manual(values = wes_palette("Royal1")) +
  scale_y_continuous(limits = c(0, 4), expand = c(0, 0)) +
  guides(fill = FALSE) +
  ggtitle("Republicans") + theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Treatment Condition", y = "")  # Add this line to set the x-axis label


rep_barplot 



rep_barplot_notitle <- means %>% filter(repdem=="Republican") %>%
  ggplot(aes(treat2, means)) +
  geom_col(aes(fill = treat2), color = "black", width = 0.85) +
  #  geom_errorbar(aes(ymin = mean_sleep - sd_sleep,
  #                    ymax = mean_sleep + sd_sleep),
  #                color = "#22292F",
  #                width = .1) +
  geom_line(data = p_value_one, 
            aes(x = x, y = y, group = 1)) +
  geom_line(data = p_value_two, 
            aes(x = x, y = y, group = 1)) +
  geom_line(data = p_value_three, 
            aes(x = x, y = y, group = 1)) +
  annotate("text", x = 1.5, y = 3.81, 
           label = p1,
           size = 3, color = "#22292F") +
  annotate("text", x = 2, y = 3.31, 
           label = p2,
           size = 3, color = "#22292F") +
  annotate("text", x = 2.5, y = 2.91, 
           label = p3,
           size = 3, color = "#22292F") +
  scale_fill_manual(values = wes_palette("Royal1")) +
  scale_y_continuous(limits = c(0, 4), expand = c(0, 0)) +
  guides(fill = FALSE) +
  ggtitle("(b)") + theme(plot.title = element_text(hjust = 0.5)) +
  labs(x = "Treatment Condition", y = "")  # Add this line to set the x-axis label


rep_barplot_notitle 
repdem_plot <- ggarrange(dem_barplot,rep_barplot)





### rerun models to separate certain vs. uncertain in the plot: 
vars_rep$certain2 <- vars_rep$certain
formula <- as.formula(paste("pool_concede_notstd ~  certain2 + control +", paste(row_names_vector, collapse = "+")))
rep2 <- lm_robust(formula, d=vars_rep)
modelsummary(rep2, stars=TRUE)

vars_dem$certain2 <- vars_dem$certain

formula <- as.formula(paste("pool_concede_notstd ~  certain2 + control +", paste(row_names_vector, collapse = "+")))
dem2 <- lm_robust(formula, d=vars_dem)
modelsummary(dem2, stars=TRUE)





demrepmodels <- list(
  "Democrats" = dem0,
      "Democrats" = dem1, 
    "Democrats" = dem2, 
  "Republicans" = rep0,
      "Republicans" = rep1,
    "Republicans" = rep2
  )



demmodels <- list(
  "Information vs. Control" = dem0,
  "Certain & Uncertain vs. Control" = dem1, 
  "Certain vs. Uncertain" = dem2
)
modelsummary(demmodels, stars=TRUE)

repmodels <- list(
  "Information vs. Control" = rep0,
  "Certain & Uncertain vs. Control" = rep1, 
  "Certain vs. Uncertain" = rep2
)
modelsummary(repmodels)

models<-list(
  "Information vs. Control" = lm_info, 
  "Certain & Uncertain vs. Control" = lm_certain, 
  "Certain vs. Uncertain" = lm_certain2
)

cmap <- c('certain2' = 'Certain \n(base Uncertain)', 'certain' = 'Certain \n(base Control)', 'uncertain' = 'Uncertain \n(base Control)', 'info' = 'Information \n(base Control)')

# Democratic models plot
demcoefplot <- modelplot(demmodels, coef_map = cmap) + labs(x = " ")

demcoefplot2 <- demcoefplot +
  scale_color_manual(values = c("black", "black", "black", "black")) +
  theme(panel.background = element_rect(fill = "grey95"),
        panel.grid.major = element_line(color = "white", size = 0.5),
        panel.grid.minor = element_line(color = "white", size = 0.25),
        legend.position = "none") +  # Remove the legend
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +  # Add a vertical line at x=0
  ggtitle("(c)") +
  theme(plot.title = element_text(hjust = 0.37)) +  # Center the title
  xlim(-0.5, 0.85)  # Set the x-axis limits
demcoefplot2

# Republican models plot
repcoefplot <- modelplot(repmodels, coef_map = cmap) + labs(x = " ")

repcoefplot2 <- repcoefplot +
  scale_color_manual(values = c("#C93312", "#C93312", "#C93312", "#C93312")) +
  theme(panel.background = element_rect(fill = "grey95"),
        panel.grid.major = element_line(color = "white", size = 0.5),
        panel.grid.minor = element_line(color = "white", size = 0.25),
        legend.position = "none") +  # Remove the legend
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +  # Add a vertical line at x=0
  ggtitle("(d)") +
  theme(plot.title = element_text(hjust = 0.37)) +  # Center the title
  xlim(-0.5, 0.85)  # Set the x-axis limits
repcoefplot2

## combine histogram with coef plots. 
dem_combine <- ggarrange(dem_barplot_notitle,demcoefplot2, ncol = 1) + ggtitle("Democrats")
dem_combine
dem_combine_with_title <- annotate_figure(dem_combine,
                                          top = text_grob("  Democrats", color = "black", face = "bold", size = 14))
dem_combine_with_title

rep_combine <- ggarrange(rep_barplot_notitle,repcoefplot2, ncol = 1)
rep_combine
rep_combine_with_title <- annotate_figure(rep_combine,
                                          top = text_grob("  Republicans", color = "black", face = "bold", size = 14))
rep_combine_with_title
repdeb_combine <- ggarrange(dem_combine_with_title,rep_combine_with_title, ncol=2)
repdeb_combine
#ggsave("repdemind_bar_coefplot.jpg",repdeb_combine)

names<- c("info" = "Information", "certain"="Certain", "uncertain"="Uncertain", "control"="Control")
#"Debt General" = "debtpreA", "Debt Taxes" = "debtpreC",
note <- paste("Non-zero covariates include in Dem. Models: ", demnotes, "Non-zero covariates include in Rep. Models: ", repnotes , collapse="")
modelsummary(demrepmodels, stars=T,coef_map = names,  notes=note )

## Tables for the Appendix (Table A3 )
#modelsummary(demrepmodels, stars=T,coef_map = names, notes=note , output="latex")
#modelsummary(demrepmodels, stars=T,notes=note , output="latex")



#### manipulation check: 











### ### ### ### ### ### 
### PETITION OUTCOME 
### ### ### ### ### ### 

sumvars <- as.data.frame(rhs)
dvsum <- vars %>% select(signpet)
treatsum <- vars %>% select(certain, uncertain, info)
sumvars <- cbind(sumvars,dvsum,treatsum)
rm(dvsum, treatsum)
rhs<- as.data.frame(rhs)
#rhs <- rhs %>%select(-partyid, -inc, -education, -RWnews, -approval_biden, -urbrur )

rhs<-as.matrix(rhs)
dv <- as.matrix(vars$signpet)

row_names_vector <- get_nonzero_coefficients(rhs, dv)

vars <- cbind(vars,rhs)

### Comparison of Info vs. No Info 
formula <- as.formula(paste("signpet ~ info + ", paste(row_names_vector, collapse = "+")))
lm_info <- lm_robust(formula, d=vars)
modelsummary(lm_info, stars=TRUE)

# Comparison of Certain and Uncertain vs. No Info 
formula <- as.formula(paste("signpet ~ certain + uncertain +", paste(row_names_vector, collapse = "+")))
lm_certain <- lm_robust(formula, d=vars)
modelsummary(lm_certain, stars=TRUE)

vars$certain2 <- vars$certain
# Comparison of Certain vs. Uncertain 
formula <- as.formula(paste("signpet ~ certain2 + control +", paste(row_names_vector, collapse = "+")))
lm_certain2 <- lm_robust(formula, d=vars)
modelsummary(lm_certain2, stars=TRUE)

models<-list(
  "Information vs. Control" = lm_info, 
  "Certain & Uncertain vs. Control" = lm_certain, 
  "Certain vs. Uncertain" = lm_certain2
)
modelsummary(models, coef_map = cmap , statistic = "{p.value} [{conf.low}, {conf.high}]")

modelplot(models, coef_map = cmap )



cmap <- c('certain2' = 'Certain (base Uncertain)', 'certain' = 'Certain (base Control)', 'uncertain' = 'Uncertain (base Control)', 'info' = 'Information (base Control)')

coefplot <- modelplot(models, coef_map = cmap) + labs(x = " ")

coefplot2 <- coefplot + 
  scale_color_manual(values = c("#C93312", "#C93312", "#C93312", "#C93312")) + 
  theme(panel.background = element_rect(fill = "grey95"),
        panel.grid.major = element_line(color = "white", size = 0.5),
        panel.grid.minor = element_line(color = "white", size = 0.25)) + 
  ggtitle("DV: Willigness to Sign a Petition") + 
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "none",  # Remove the legend
        legend.background = element_rect(fill = "white", colour = "black"), 
        legend.key = element_rect(fill = "white", colour = "white")) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black")

coefplot2
#ggsave("pool_petition_coefplot.jpg",coefplot2)









